function y = DistP2Segment(P,Seg)

% required that Seg includes 2 distinct points
% this version is SLOW in calculation, computational expensive
% refer to calculation in note p16

a = (Seg(2,2)-Seg(2,1))/(Seg(1,2)-Seg(1,1));   
b = (Seg(1,2)*Seg(2,1)-Seg(1,1)*Seg(2,2))/(Seg(1,2)-Seg(1,1));

dBC = [Seg(1,2)-Seg(1,1);Seg(2,2)-Seg(2,1)];

a2 = dBC(1);
b2 = dBC(2);
c2 = P(1)*a2+P(2)*b2;

A = [a2 b2; a -1];
det_A = A(1,1)*A(2,2)-A(2,1)*A(1,2);
B = [c2; -b];      % B here not point B but is B matrix in Cramer's rule
A1 = [B A(:,2)];
det_A1 = A1(1,1)*A1(2,2)-A1(2,1)*A1(1,2);
A2 = [A(:,1) B];
det_A2 = A2(1,1)*A2(2,2)-A2(2,1)*A2(1,2);

H = [det_A1/det_A; det_A2/det_A];

minx = min(Seg(1,:));
miny = min(Seg(2,:));
maxx = max(Seg(1,:));
maxy = max(Seg(2,:));

if (minx<=H(1))&&(H(1)<=maxx)&&(miny<=H(2))&&(H(2)<=maxy)
    y = norm([H(1)-P(1) H(2)-P(2)]);
else 
    dP2B = norm([Seg(1,1)-P(1) Seg(2,1)-P(2)]);
    dP2C = norm([Seg(1,2)-P(1) Seg(2,2)-P(2)]);
    y = min(dP2B, dP2C);
end

end